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Abstract 



A subroutine for very-high-precision numerical solution of a class of ordinary differential equations is provided. For given evalua- 
tion point and equation parameters the memory requirement scales linearly with precision P, and the number of algebraic operations 
scales roughly linearly with P when P becomes sufficiently large. We discuss results from extensive tests of the code, and how one 
for a given evaluation point and equation parameters may estimate precision loss and computing time in advance. 
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q PROGRAM SUMMARY 

y—{ Manuscript Title: High precision series solution of differential 
equations: Ordinary and regular singular point of second order ODEs. 
' 1 ' Authors: Amna Noreen, Kare Olaussen 
^2 Program Title: seriesSolveOdel 
I Journal Reference: 
, Catalogue identifier: 

Licensing provisions: none 

S Programming language: C++ 
Computer: PC's or higher performance computers 
Operating system: Linux and MacOS 
_^ RAM: Few to many megabytes (problem dependent) 

^ Number of processors used: 1 
\Q Keywords: Second order ODEs, Regular singular points, Ordinary 
O^l points, Frobenius method. 

Classification: 2.7 Wave functions and integrals, 4.3 Differential 
equations. 

External routines/libraries: CLN - Class Library for Numbers (TJ 
O built with the GNU MP library 0, and GSL - GNU Scientific 
O^l Library (3) (only for time measurements). 

Subprograms used: The code of the main algorithm is in 
^ the file seriesSolveQdel.ee, which #include the file 
•rH checkForBreakOdel.ee. These routines, and programs using 

them, must #include the file seriesSolveOdel . cc. 



executing the recursion 



^ Nature of problem: The differential equation 



tf_ 

dz 2 



■V-d v+ v_ 



dz 



1 N 

<Mz) + - £ v„ z " <Kz) = 0, (1) 



is solved numerically to very high precision. The evaluation point z 
and some or all of the equation parameters may be complex numbers; 
some or all of them may be represented exactly in terms of rational 
numbers. 
Solution method: 

The solution i//(z), and optionally tft'(z), is evaluated at the point z by 
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A m+ ,(Z) = ~. : : r V V„(z) A m _„(z), (2) 

(m + I + v - v + )(m + 1 + v - y_) 

(3) 



^ m+l \z) = ^ m, {z)+A m+l {z), 



to sufficiently large m. Here v is either v + or v_, and V„(z) = v„ z" +1 . 
The recursion is initialized by 



A- n (z) = 6„ a z v , torn = 0, 1,...,N 
<A (0) (z) = A (z). 



(4) 
(5) 



Restrictions: No solution is computed if z = 0, or s = 0, or if v = 
y_ (assuming Rev + > Rev_) with y + - y_ an integer, except when 
v + - v_ = 1 and vo = (i.e. when z is an ordinary point for z~ y ~ (Hz))- 
Running time: On an few years old Linux PC, evaluating the ground 
state wavefunction of the anharmonic oscillator (with the eigenvalue 
known in advance), cf. equation 0, at y = VlO to P = 200 decimal 
digits accuracy takes about 2 milliseconds, to P = 100000 decimal 
digits accuracy takes about 40 minutes. 
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1. Introduction 

Modelling and analysis of many problems in science and en- 
gineering involves the solution of ordinary differential equa- 
tions, sometimes in the domain of complex numbers. For prac- 
tical use such solutions must usually be computed numerically 
at some stage. In some cases it may be necessary, useful, or 
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interesting to do this to much higher precision than provided by 
standard equation solvers (or routines for evaluating standard 
functions). 

Some examples of cases where numerical calculations have 
been used to inspire or check analytic conjectures and proofs 
are the works by Bender and Wu G] and Zinn-Justin and 
Jentschura [2]. With access to very accurate numerical results 
the opportunities for such explorations increases. Access to 
essentially exact results are also useful for analyzing the be- 
haviour of approximation schemes, as in the work by Bender et 
al G) and more recently by Mushtaq et al [4|. 

We have implemented and investigated the algorithm ([2] [3) 
for solving equation ([T} to very high precision, and believe the 
C++ function seriesSolveOdel may be of use or interest to 
others. An early version of this code has been used to solve 
eigenvalue problems like the anharmonic oscillator and the dou- 
ble well potential, 
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2 ■ y 



i//(y) = e„ if/(y), 

2 



dy 2 



(i-r)' 



(6) 



(7) 



to very high precision. In reference Q the ground state eigen- 
energy eo of (|6} was found to 1 + million decimal digits pre- 
cision, the excited state £50000 was solved to 50000 + decimals, 
and the lowest even, eq+, and odd parity, eo-, eigenvalues of Q, 
with s = 1/50000, was found to 30000 + decimals. Equations 
(|6] |7ji are transformed to the form ([T]i by introducing z — y 2 , 
leading to v_ = and v + = | and V2 = 7. This further gives 
Vo = -\s„ for equation (^J, and Vo = |(1 - e ncr ), vi = -\ for 
equation Q. I.e., the eigenvalue parameter enters in the coeffi- 
cient Ao(z) of equation (|8j. 

The eigenvalue condition for these problems is that the wave 
function should vanish as y — > ±00, a condition which cannot be 
imposed numerically. However, an asymptotic analysis of the 
behaviour of the wavefunction as y — > +00 allows us to replace 
it with an equivalent Robin boundary condition at finite y. The 
latter cannot be computed exactly, but to sufficient accuracy for 
any desired precision. In fact, if we make y large enough it 
suffices to use a Diriclet boundary condition. 

In reference (6] it was demonstrated that the wavefunction 
normalization integrals can be computed to comparable preci- 
sion, again using an early version of our code. 

In the rest of this paper we provide examples of how this 
code can be used, and some analysis of its behaviour. We do 
not focus on specific areas of applications, but would like to 
mention that very-high-precision computations of Green func- 
tions and functional determinants are possible applications. The 
code can evaluate ifr(z) for complex values of z, and allow for 
complex parameters in the differential equation. 



2. Basic use 



The function seriesSolveOdel is declared as 



Function declaration 

bool seriesSolveOdel (OdeResultsfe results, const 
OdeParamsfe params) 

The function returns true if the calculation completed nor- 
mally and false if the calculation was aborted. Function argu- 
ments and options are collected in a structure OdeParams, with 
the function value i/Kz), optionally i//(z), and various diagnostic 
results returned in a structure OdeResults. The definitions of 
these structures are listed in the Appendix at the end of this pa- 
per. A code snippet illustrating the use of seriesSolveOdel 
is 

Example use 

II Better start from default parameters 
OdeParams params = def aultOdeParams ; 
// Change parameters as needed 
params. prec = f loat_f ormat (500) ; 
params. z = complex (27/2 , 43/7); 
params. dAlso = true; 
OdeResults results; 

// NB! Must allocate space for tft(z) a if/'(z) 
cl_N fu[2] ; 
results. fu = fu; 

if (seriesSolveOdel (results , params) ) { 
cout << results . fu [0] << endl; 
cout << results . fu [1] << endl; 

} 



3. Computational accuracy 

Our solution is found by a brute force summation, 



M 



M 



a m z 



2 a m z v+m = Y, AJz) - 



(8) 



Since equation ([T} has no singularities in the finite z-plane, ex- 
cept perhaps z — 0, the sum in guaranteed to have an infinite 
radius of convergence. However, intermediate terms in the sum 
may be very large although the final result is small; hence there 
may be huge cancellations, leading to significant loss of accu- 
racy when z is large. 

The computations are performed with high-precision float- 
ing point numbers, with precision regulated by the parameter 
params .prec as shown in the code snippet above. The value 
given is the intended precision P in decimal digits, but since 
memory for floating point numbers is allocated internally (in 
CLN at high precision) in chunks of 64 bits the actual preci- 
sion is usually somewhat higher — increasing in steps of about 
19 * 64 • log 2/ log 10 decimal digits. 

Denote the actual precision used in computation by M. 
I.e., each non-zero floating point number x is represented as 
(-l) J 2 m /, with the mantissa / (1 < f < 1) given to a precision 
of M bits. This means that the potential roundoff error in x is 
2»n-M-i j e> tne laj-ggst t erm A m (z) in the sum §8§ has the rep- 
resentation (-1) J 2 A / it may contribute a roundoff error 2 A ~ M ~ l 
to Due to the recursion relation |2]l roundoff errors may be 



further amplified (or partially cancelled), but it is a reasonable 
hypothesis to use A for an estimate of the evaluation error. 




Figure 1: The real evaluation error, s r = \ipp(z) - caused 
by roundoff is strongly correlated with the largest term A max in 
the series expansion. Here if/p(z) is the value obtained when 
evaluating the series to the intended precision of P decimal dig- 
its, while <Jf E (z) is the "exact" value (here actually the value 
obtained for P = 1040). 

We have tested this hypothesis by running a large number 
(200000) of evaluations with randomly chosen parameters, and 
investigated the correlation between A and the real evaluation 
error s r = \iffp(z) - ^e(z)\ for various precisions P. Here i// E (z), 
representing the exact value, is found by doing the computation 
with a precison P E reasonably larger than all the others. As can 
been seen in figure [I] the correlation between A and s r is good 
compared to the accuracies in question. 

For diagnostic purposes the value of A is returned by 
seriesSolveOdel in the variable results .maxAExponent. 
The corresponding value for iff'(z) is returned in the vari- 
able results . maxAdExponent (the values of m where the 
maxima occur are also returned). Based on these values 
and the actual precision M the estimated errors in deci- 
mal digits are returned in the variables results . lgErrorF 
and results . lgErrorFd. Their exact values are (A - 
M)log2/log 10 + G for iff, and (A' - M)log2/log 10 + G' for 
iff', where the numbers G = 4.30 and G' = 3.02 are empirically 
choosen "guard digits" to avoid underestimating the error (too 
often). 

These estimates are accurate to a handful of decimal digits 
as shown by the A = lg |s r | - lg \e e \ histogram in figure [2] We 
have found such histograms to be independent of computational 
precision P, and (with the choosen value of G - G') also the 
same for tf/(z) and if/'(z). The histogram is taken over 200000 
evaluations with N choosen randomly between 1 and 4, the real 
and imaginary parts of s randomly from the set {— 1, — 4, i, 1}, 
the real and imaginary parts of v ± randomly between -10 and 
10, the real and imaginary parts of v„ randomly between -5 and 



5, and the real and imaginary parts of z randomly between -20 
and 20. 



Histograms of error differences for tf>(z) 
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Figure 2: The curve for A = lg |e r | - lg \e e \ shows a histogram of 
differences between the real evaluation error s n and an estimate 
based on the largest term A max in the sum (fSJ. The histogram 
does not depend on computational precision P, but individual 
differences (with fixed z and equation parameters) varies with 
P in a fluctuating manner. The curve for A - A shows how 
differences computed at precision P — 20 (A) correlates with 
those computed at P = 200, 500, 1000. 

The difference A is caused by an essentially unpredictable 
roundoff error, amplified by a recursion relation which depends 
on z and parameters of the differential equation. As can be seen 
from figure [2] the ratio between the real and estimated errors 
varies between almost 10 5 and 10~ 8 . Although this variation is 
large it is still less than the discrete steps by which the actual 
precision is increased. 

We have investigated how the difference A computed at dif- 
ferent precisions (but for the same evaluation point z and equa- 
tion parameters) are correlated. This is shown in the A - A 
histogram in figure [2J where A refer to a computation with in- 
tended precision P = 20, and A to computations with P = 200, 
500, and 1 000 (where each P gives the same looking his- 
togram). As can be seen the correlations are stronger than for a 
single A, but there are still wide tails. 

In conclusion, the largest term A max in the series (|8j, or the 
associated integer A, provides a good estimate of the evalua- 
tion error, but the real error may still differ by several orders 
of magnitude. As somewhat better empirical estimate can be 
obtained by first computing the real error at low P (where it is 
computationally inexpensive) and assuming 

lg se = lg £ r - (M - M) log 2/ log 10 + 2. (9) 

Here s r is the real error at M bits of actual precision, with se 
the estimated error at M bits of actual precision. 
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4. Comparison with exactly known Wronski determinant 



In the previous section we assumed that seriesSolveOdel 
would compute accurate results at large intended precision P, 
but this was not really verified. One check is to compare its 
results with the exactly known Wronski determinant, 



W(z) = ifr Vt (z)KJz) ~ <Av_(z)<A; + (z) = (v- - v+)z 



v + +v_-l 



(10) 
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Fraction 



• *• 



VgSW r 



530 



-525 



-520 



510 



-505 -500 



Figure 3: The real errors 8W r = |W exact (z) - W num (z)\ when 
aiming to compute the Wronski determinant W(z) to an abso- 
lute accuracy of 10~ 5l)0 . The histogram is taken over 150000 + 
random evaluation points z and equation parameters. 



The numerically computed determinant is estimated to have 
an error of magnitude 

SW e = max (|^ v+ | &f/ y _, \t/r y _\ 6if/ Y+ , \\f/ v \d^,W^\6i{i y ], (11) 

where f.i. Sif/ V is the estimated magnitude of error in i// v (all 
quantities evaluated at z). An estimate of the loss of precision 
can be made by a calculation at low(er) P, and used to choose 
the appropriate value of params .prec for a desired final pre- 
cision in W(z). Figure [3] shows a histogram of how this works, 
tested on a large number of random evaluation points and equa- 
tion parameters. In most cases the real precision is reasonbly 
close (always better) than the one aimed for, but sometimes it 
turns out to be much better. This may occur when the low pre- 
cision calculation overestimates the magnitude of tf/(z) or tfr'(z). 
However, as shown in figure [4] the real error 5W r in the numer- 
ically computed determinant is always reasonably close to the 



final estimate 6W e based on equation (Hi with all quantities 
computed at high precision. 
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Figure 4: Histogram of differences between the real error 6W r 
in the numerically computed Wronski determinant, and the es- 
timated error 6W e based on ( flT) . 

The conclusion is that the calculated Wronski determinants 
are correct within the expected accuracies, at least for compu- 
tations at sufficiently high precision (how high may depend on 
the evaluation point z and equation parameters). 

5. A priori accuracy estimates 

Although seriesSolveOdel monitors the largest term in 
the sum <|8j to estimate the accuracy of the computed results, 
it is desirable to predict the behaviour of the sum in advance. 
One way to do so is by analysing the behaviour of the recur- 
sion relation (|2j. Quite detailed and interesting results can be 
found in simple cases, but the analysis becomes unmanageable 
in general. We instead make the hypothesis that the terms in 
the sum <JsJ for large z = xe llf is strongly peaked (in absolute 
value) around some m = m, and that there exist values of ip for 
which there are little cancellation between the large terms. I.e., 
we assume that 



max (fr(ie 1( °) » 



(12) 



for positive x. We may use the WKB-approximation to estimate 
the left hand side. For analytic treatment we first neglect the 
slowly varying algebraic prefactor of the WKB-approximation, 
and a similar correction to the relation ( |12| l. Such corrections 
can be included in a numerical implementations. 
Define, for positive u, 



S(u) = max log |i^(e 



max Re 

<p 



f 

Jo 



Q(t)dt 



(13) 



where Q(t) is found from the differential equation (JTji. We then 
have the relation 



S (u) = log (\dm\) + (v + m) u, 
d 



dm 



log(|fl m |) 



(14) 
(15) 
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The last equation follows from the maximum condition. We 
recognize (14 15 i as a Legendre transform [|8'|,|5], IfTOll l. By 



inverting this transform we find 



Note that (19 20 1 give, in parametric form, an estimate of all 



coefficients \a m \, not only those corresponding to a maximum 
value. For c — an explicit representation is easily found to be 



m-—S(u), log (|flm|) = S(u) - (v + to) u, (16) 
du 

which provides an a priori order-of-magnitude estimate of the 
coefficients a m , and hence of (i) the accuracy loss due to numer- 
ical roundoff, and (ii) the number of terms M required in <[8j for 
a desired final precision. 



5.1. Example 1: Anharmonic oscillators 
Consider the equation 

-^(y) + (/ + c 2 ) 2 ¥(y) = 0. (17) 
For large y the typical solution behaves like 

«F(y) ~ e s y3+cZy , (18) 



neglecting the slowly varying prefactor. Equation ( 17 1 can be 
transformed to the form ([T} by introducing x = y 2 , *F(y) = ifr(x). 
Hence, with x — y 2 = e" 



S(u)= i(e5" + 3c 2 e5"), 



which gives 



m = i(e5" + c 2 e5"), (19) 
log(|o*|) = (§- \u)& u +c 2 (l - I«)ei". (20) 

Numerical and estimated coefficients 
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Figure 5: Comparison of numerical coefficients a m (points) 
with estimates (full-drawn lines) based on ( p~9| [20] > and ( [28] [29) . 
The estimates of log \a m \ are accurate up to corrections which 
depend logarithmically on to. 



log \a m \= -m (1 - log 2m) . 



(21) 



This is plotted as the lower curve in figure [5] It fits satisfactory 
with the high-precision coefficients generated numerically, but 
there remains a correction which depends logarithmically on to. 
For nonzero c the parametric representation provides equally 
good results, as shown by the upper curve in figure [5] 

The conclusion of this example is that we expect the largest 
term of the power series to be 



max|A m (jt)| ~ e5 



(22) 



neglecting a slowly varying prefactor. Further, the maximum 
should occur at 



\ (x 312 + c 2 x 



l/2\ 



(23) 



Finally, estimates like equation (21 1 for the coefficients a„, may 
be used to predict how many terms Ad we must sum to evaluate 
iff(x) to a given precision P, based on the stopping criterium 



IgmIx^ < 10" 



(24) 



As can be seen in figure [6] the agreement with the actual num- 
ber of terms used by seriesSolveOdel is good, for case of 
equation ( [T7] i with c — 0, in particular for high precision P. But 
note that a logarithmic scale makes it easier for a comparison to 
look good. 



Predicted and calculated summation lengths 




Figure 6: This figure compares the a priori prediction, based 



on equation (21 1, of the number of terms M which must be 
summed in order to evaluate for c = to a desired 
precision P with the actual number of terms computed by 
seriesSolveOdel. 
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5.2. Example 2: Double well oscillators 



5.3. Logarithmic corrections 



Next consider the equation 



Ratio of numerical to estimated coefficients 



dy 



- x ¥(y) + (y 2 -c 2 fw(y) = 0. 



(25) 



For large y the typical solution behaves like 



(26) 



neglecting the slowly varying prefactor. Equation ( 25 i can be 
transformed to the form ([!} by introducing x = y 2 , ( F(y) = t//(x). 
Hence, with x - y 2 = e" 

S(u) = max iRe(e5 ( " +i ^ - 3c 2 eJ ( " +i) ^) . 



The maximum occurs for cos ^tp = -j (l + c 2 e "J 
e" > |c 2 , and for cos lip = -1 otherwise. This gives 



c 2 e" /2 - ±e 3 " /2 fore" < \c 2 , 
S< " } ~ ' ±(e" + c 2 ) 3/2 fore" > \c 2 . 



1/2 



when 



(27) 



This implies that 



(28) 



log (la^l) 



±e M/2 (c 2 -e") fore"<|c 2 , 

ie" (e" + c 2 )' /2 fore">|c 2 , 

(l - ±h)cV /2 - (i - i M )e 3 " /2 for e" < ±c 2 , 

[(i - it<)e" + ±c 2 ] (e" + c 2 )' /2 for e" > ±c 2 . 

(29) 



This representation compares fairly well with the numerically 
generated coefficients, as shown by the middle curve in figure[5] 
However, in this case the coefficients a m have a local oscillating 



behaviour. The representation (28 29 1 should be interpreted as 
the local amplitude of this oscillation. 

The conclusion of this example is that we expect the largest 
term of the power series to be term of the series to be 




Figure 7: The ratio between the numerically gene rated \a m \ and 
their estimated values a„ based on equation ( 16 1 with S(u) re- 
placed by S e ff(u). The important feature here is that the ratios 
are essentially constant, not which constant, due to an overall 
normalization constant which we have not attempted to esti- 
mate here. 

The agreement between the a priori estimated magnitude a (e) 
and the numerically generated values \a m \ looks quite good in 
figure [5] This is partly due to the logarithmic scale; on closer 
examination the coefficients are seen to differ by many orders 
of magnitude. The agreement can be significantly improved 
by (i) taking into account the prefactor Q(uY 1 ^ 2 in the WKB- 
expression ( p~3] > for the left hand side of \\2) , and (ii) changing 
the right hand side of ( 12 1 as 



la- \x v+m = e I (" 1 )+0 /+m )" 



V27r5"(M)e" 



(m)+(v+m)u 



The latter replacement takes into account that the main contri- 
bution to the sum over m comes from a range of values around 
m, approximates this contribution by a gaussian integral, and 
uses the fact that s"(/«)~' = -S"(u), These two improvements 
amounts to a change S(u) — > S e ff(u). 



Implementing these changes for the c - case of ( 17 1, taking 
into account that only coefficients 03,, + 0, the estimate pi) can 
be improved to 



max |A m (x)| 



(30) 



neglecting the slowly varying prefactor. Further, the maximum 
should occur at 



; (x + c 2 ) 



1/2 



\x" 2 + \c 2 x 1 ' 2 . 



(31) 



log a! 



(e) 



I (m + I) [l - log(2m + §)] - i log f . (32) 



Similar, algebraically more complicated, improvements can be 
made for c > 0. As shown in figure [7] the improved estimates 
compare very well with the numerically generated coefficients, 
but in case of locally oscillating a m the (smooth) estimate a$ 
should be interpreted as the local oscillation amplitude (as il- 
lustrated by the (y 2 - 5 2 ) 2 -case in figure]^. 



6. Concluding remarks 

For equations without singular points in the finite plane our 
code can be used for expansion around any point £o in the com- 
plex plane. To shift from one expansion point to another one 
just has to rewrite the parameters v„, and let z denote the dis- 
tance from £q. This allows for analytic continuation of the solu- 
tion, which becomes quite easy since the full solution is deter- 
mined by just the two parameters tff((o), ^'(£0) (in addition to 
the differential equation). 

The strategy of using a sequence of series expansions, each 
with a small parameter z, has been used by Haftel etal Bill . The 
advantage is that each summation requires fewer terms in the 
series, and may lead to less loss of precision caused by roundoff 
errors. The cost is of course that one has to do several sums, and 
one may also loose symmetries like the y — > —y symmetry in 
equations |6] |7J. The latter leads to more algebraic operations 
per recursion step. 

The optimal strategy may depend on the problem. If we are 
only solving equation (|6]l for the ground state eigenvalue £0, 
this is basically determined by the condition that the asymptotic 
behaviour of the solution switches very rapidly between e v ^ 



and 



This behaviour is not affected much by roundoff 



errors. Consider the question is whether it is faster to evaluate 
e v / 3 by a single series expansion, or by k expansions with yk = 



y/k. By combining equations (21 24 1 one finds that each sum 
requires about Mk terms for a given precision P, where Af* 
satisfies the equation 

|Afc (1 - log 2M k ) + 2M k logfy/k) * -P log 10, (33) 

which is best solved numerically. Consider f.i. the case of 
y = Vl78 and P - 10 5 . As can be seen from figure ^ about 
M = M\ = 10 5 terms have to be summed to obtain the desired 
precision. With k — 2 only about M2 = 67 500 terms has to 
be summed, but since this has to be done twice the total effort 
becomes larger. The situation is similar for other values of y 
and P. 

In other cases, like highly excited states of |6]l or all states of 
0, there is a loss of precision due to roundoff. This changes 
how the number of terms M and the actual precision M vary 
with y. The latter is most important since the time per multi- 
plication increases somewhat faster than quadratic with M. In 
such cases a sequence of analytically continued evaluations are 
clearly advantageous; optimization of the number and size of 
steps requires some prior knowledge of the coefficients a„ and 
the behaviour of the solution lfT3l . If one needs to evaluate 
the solution at a sequence of points, as when calculating the 
normalization integral |6), analytic continuation would also be 
preferrable. 

The routine seriesSolveOdel does not allow for analytic 
continuation in the presence of a regular singular point, since 
the transformed equation belong to a different class. We have 
developed and are testing code for a more general class of equa- 
tions (as hinted by our naming scheme), which we intend to 
submit real soon. This code allow for translations (or more gen- 
erally Mobius transformations) to a new expansion point (q. It 
can f.i. be used to solve Mathieu and Mathieu-like equations. 



We have made extensive tests of the submitted code, which 
appears to be robust and perform according to theoretical ex- 
pectations. As illustrated, surprisingly accurate a priori esti- 
mates of the series to be summed can be made by using the 
WKB approximation in combination with Legendre transforms. 
This is useful for estimating precision and time requirements in 
advance. In the general case the WKB integral and Legendre 
transformation must be computed numerically. We have devel- 
oped and tested code for this purpose lfT2l[T3l . 
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Appendix 



Definition of the OdeParams structure 

struct OdeParams { 

// Parameters defining the ODE 

z; // Evaluation point 

nuP; // v+ 

nuM; // v_ 



cl_N 
cl_N 
cl_N 
cl_N 
int 
cl_N* 



orderN; 



bool nuPlus ; 

bool dAlso; 

cl_I emmMax; 

cl_I emmTooLarge ; 



// Order of polynomial potential 

// Pointer to array v[orderN+l] 
// Options regulating the computation 

f loat_f ormat_t prec; // Wanted computational precision 

// If true, compute v = v + solution 

// If true, compute derivative also 

// Stop summation when emm=emmMax 

// Stop if emm >= emmTooLarge 
// Options to print double precision coeff info to stdout 

bool writeParams; // Write to (temporary) file? 

// Delete above file at normal exit? 

// Option to print all log(abs(A„,)) 

// Option to print log(abs(Re(A„,))) 

// Option to print log(abs(Im(A„,))) 

// Option to print arg(A,„) 

// Option to print Re(A,„) 

// Option to print Im(A„,) 

// Print all log(abs((v + m)A m /z)) 

II Print log(abs(Re((v + m)A m )/z)) 

II Print log(abs(Im((v + m)A m /z))) 

II To print arg((v + m)A m )/z 

II To print Re((v + m)A m /z) 

II To print Im((v + m)A m /z) 



deleteParams ; 

printLogAbsA ; 

printLogReA; 

printLoglmA; 

printArgA ; 

printReA; 

printlmA; 

printLogAbsAd; 

printLogReAd; 

printLoglmAd; 

printArgAd; 

printReAd; 

printlmAd; 



bool 
bool 
bool 
bool 
bool 
bool 
bool 
bool 
bool 
bool 
bool 
bool 
bool 

// Options to print full precision coefficients to stdout 

bool coutA; // Option to print A m 

bool coutReA; // Option to print Re(A,„) 

bool coutlmA; // Option to print Im(A,„) 

bool coutAd; // Option to print (v + m)A,„/z 

bool coutReAd; // Option to print Re((v + m)A m /z) 

bool coutlmAd; // Option to print Im((v + m)A m /z) 



struct OdeResults 
cl_N* fu; 
int maxAExponent ; 
int maxAdExponent 
int emmAtAMax; 
int emmAtAdMax; 
int lengthOfSum; 
double timeUsed; 
double lgErrorF; 
double lgErrorFd; 
int returnStatus ; 



}; 



Definition of the OdeResults structure 
{ 

// Pointer to pre-created array fu[2] 

// Largest term in function series 

// Largest term in derivative series 

// m where abs(A m ) is largest 

// m where abs((v + m)A m ) is largest 

// Number of terms summed 

// Time used for evaluation 

// Estimated error in function 

// Estimated error in derivative 



// ■ 
// ■ 
// ■ 
// ■ 
// 1: 
// 2: 



s = 
; = 

(vp - vm) is integer ! = 1 
emm >= emmTooLarge was reached 
All abs(A,„) < Estimated error 
emm = emmMax was reached 
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